function [Grid,TrMat]=Rouwenhorst(N,Rho,Mean,Sig,Skew)

%--------------------------------------------------------------------------
% Note:
%   Discrete Grid: [m-d,m-d+1*2d/(N-1),m-d+2*2d/(N-1),...,m+d]
%   E[y_n]          =   m+d*( (q-p)/(2-p-q) )
%   Var[y_n]        =   4*d^2/(N-1)*( (1-p)(1-q)/(2-p-q)^2 )
%   Corr[y_n,y_n+1] =   p+q-1
%   Skew[y_n]       =   (p-q)/sqrt( (1-p)(1-q) )/sqrt(N-1)
%   Kurt[y_n]       =   ( -2+(p-q)^2/( (1-p)(1-q) ) )/(N-1)
%% Preliminary
% CHECK IF abs(rho)<1 
if abs(Rho)>=1
    error('The persistence parameter Rho is not smaller than 1!');
end

% CHECK IF n_R IS AN INTEGER GREATER THAN ONE.
if N<1.50001
    error('The number of grid points N must be greater than one.');  
end

% CHECK IF n_R IS AN INTEGER.
if mod(N,1)~=0 
    warning('The number of grid points is not an integer! Rounded value will be used!')
    N=round(N);
end

%% Construct the Grid and Transition Matrix
%--------------------------------------------------------------------------
% Solve q
%--------------------------------------------------------------------------
% Symmetric Process
if Skew==0
    q       =   (1+Rho)/2;
else
    Target  =   sqrt(N-1)*Skew;
    TempFun =   @(q)(1+Rho-2*q)/sqrt((1-q)*(q-Rho))-Target;
    q       =   VecFun_Bisection(TempFun,Rho+eps*10,1-eps*10,1e-6);
end
p       =   1+Rho-q;
d       =   Sig*sqrt( (N-1)*( (2-p-q)^2/(1-p)/(1-q) )/4 );
m       =   Mean-d*(q-p)/(2-p-q);

Grid    =   m+d*linspace(-1,1,N)';
TrMat   =   [p 1-p;1-q q];
ii      =   2;
while ii<N
    Mat_1 = [TrMat,zeros(ii,1);zeros(1,ii+1)];
    Mat_2 = [zeros(ii,1),TrMat;zeros(1,ii+1)];
    Mat_3 = [zeros(1,ii+1);TrMat,zeros(ii,1)];
    Mat_4 = [zeros(1,ii+1);zeros(ii,1),TrMat];
    TrMat = p*Mat_1+(1-p)*Mat_2+(1-q)*Mat_3+q*Mat_4;
    TrMat(2:end-1,:) = TrMat(2:end-1,:)/2;
    ii = ii+1;
end
TrMat = TrMat';
TrMat = bsxfun(@rdivide,TrMat,sum(TrMat));